library(pROC)
library(PRROC)
library(caTools)
###########################################################
###GLMs: With Interaction (Figure 5, A6-A7, and Table 3)###
###########################################################

#predicted probabilities
scores<-pred.glm08$fit

#true binary outcomes
binary.labels <- ss.dat.f08$civconf==1

#threshold
threshold<-.5

#classification stats
tp <- sum( (scores > threshold) & binary.labels )
sensitivity <- tp / sum(binary.labels)
tn <- sum( (scores <= threshold) & (! binary.labels))
specificity <- tn / sum(!binary.labels)
fp<-sum( (scores > threshold) & (! binary.labels))
fn<-sum( (scores <= threshold) & binary.labels)
False.Positive.Rate<-1-specificity
False.Negative.Rate<-fn/(fn+tp)
PPV<-tp/(tp+fp)
F1<-2*((PPV*sensitivity)/(PPV+sensitivity))
F2<-(5*(PPV*sensitivity))/(4*(PPV+sensitivity))
Accuracy<-(tp+tn)/(tp+tn+fp+fn)

#print results
sensitivity
specificity
False.Positive.Rate
False.Negative.Rate
Accuracy

#Calculate AUC and Area under PR curve
colAUC(pred.glm08$fit,ss.dat.f08$civconf, plotROC=FALSE) 
pr.glm.f <- pr.curve(scores.class0 =pred.glm08$fit, weights.class0 =ss.dat.f08$civconf, curve=T)
#Plot PRC for Figure 5
jpeg("glmprc2.jpeg", width = 6, height = 6, units = 'in', res = 500)
plot(pr.glm.f)
dev.off()
#PRC AUC
pr.glm.f[[2]]
#Store all information 
Store.2.glm <-cbind(pr.curve(scores.class0 =pred.glm08$fit, weights.class0 =ss.dat.f08$civconf)[[2]],colAUC(pred.glm08$fit,ss.dat.f08$civconf, plotROC=FALSE),PPV,sensitivity,F1,F2,Accuracy)
#Table 3
Store.2.glm

###Plot ROCs (Figure A6 and A7)
#No interaction (Figure A6)
jpeg("glmrocs.jpeg", width = 6, height = 6, units = 'in', res = 500)
roc.glm.f08s <- plot.roc(ss.dat.f08$civconf, pred.glm08s$fit,
                         main="GLM (No Maize and Interaction)",  xlab="False positives",ylab="True positives",
                         print.auc=TRUE, print.auc.x=0.5, print.auc.y=0.1, ci=TRUE)
dev.off()

#With interaction (Figure A7)
jpeg("glmroc.jpeg", width = 6, height = 6, units = 'in', res = 500)
roc.glm.f08 <- plot.roc(ss.dat.f08$civconf, pred.glm08$fit,
                        main="GLM",  xlab="False positives",ylab="True positives",
                        print.auc=TRUE, print.auc.x=0.5, print.auc.y=0.1, ci=TRUE)
dev.off()